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ABSTRACT 

The properties of a preconditioned, coupled, strongly im- 
plicit finite-difference scheme for solving the compressible 
Navier-Stokes equations in primitive variables are investi- 
gated for two unsteady flows at low speeds, namely the 
impulsively started driven cavity and the startup of pipe flow. 
For the shear-driven cavity flow, the computational effort was 
observed to be nearly independent of Mach number, espe- 
cially at the low end of the range considered. This Mach 
number independence was also observed for steady pipe 
flow c aintiaiinns; however, rather different conclusions were 
drawn for the unsteady calculations. In the pressure-driven 
pipe startup problem, the compressibility of the fluid began 
to significantly influence the physics of the flow development 
at quite low Mach numbers. The present scheme was ob- 
served to produce the expected characteristics of completely 
incompressible flow when the Mach number was set at very 
low values. Good agreement with incompressible results 
available in the literature was observed. 


INTRODUCTION 

Over the past few years a more complete understand- 
ing of the behavior of numerical algorithms for solving the 
time-dependent, compressible Navier-S tokes equations at low 
Mach numbers has evolved. Until this understanding devel- 
oped, most algorithms designed for compressible flows were 
observed to become very inefficient, or inaccurate, or both 
at low Mach numbers. This is unnecessary because it is 
well known that the physics itself is no more complex just 
because the Mach number is low. That is, for most flows, 
no important changes would be observed if the Mach number 
were reduced from, say, 02. to 0.01, if all other dimensionless 
parameters of the flow remained the same. The difficulty 
must be due to an inappropriate structure of the algorithm. 
This problem and proposed remedies have been addressed by 
several investigators including Turkel 1 - 2 , Feng andMerkle 3 , 
Peyret and Viviand 4 , and Choi and Merlde. 5 

Previous investigators have adopted different points of view 
but have generally agreed that the hyperbolic time-dependent 
Navier-Stokes system becomes “stiff” at low Mach numbers 
of greatly different signal speeds (convective and 

acoustic) which result in large differences in the magnitudes 
of the eigenvalues of the system. Explicit schemes especially 


are restricted to using very small time steps as dictated by 
stability based on acoustic speeds. Fluid particles, on the 
other hand, may require several of these small time steps to 
traverse a computational cell. As pointed out by Choi and 
Merkle 5 difficulties allow Mach numbers with implicit fac- 
tored schemes may be due in large part to factorization errors 
unless the time step is maintained quite small. The alterations 
made to the numerical formulation to overcome the stiffness 
problems of the hyperbolic system has become known as 
“preconditioning” which will be discussed in physical terms 
in a section to follow. Despite this name which may be 
adequately descriptive from some points of view, the precon- 
ditioned system is merely one with a “modified” time term. 
This can be made to appear in the equations as a new matrix 
multiplying the time term in the vector form of the system 
of equations. To solve steady problems, the preconditioning 
can be accomplished by altering the time derivative term in 
the equations. To solve unsteady problems, the physical time 
term can be left intact and an additional “pseudo-time” term 
of a particular form added to the equations. These altered or 
added time terms change the nature of the hyperbolic problem 
which is being advanced in a “new” or “pseudo” time variable. 
Because all the altered time terms vanish at c onverg ence (to 
steady state or at each physical time step), no approximation 
is bang injected into the governing equations. Essentially, in 
the low Mach number limit, the preconditioning enables the 
numerical model to smoothly bridge the gap between a fully 
compressible formulation and one' that can, for isothermal 
flows, be viewed as a pseudo-compressible formulation 6 of 
the incompressible equations. 

Low speed flows can, of course, be computed by using a 
completely incompressible formulation. One of the objectives 
of the present paper is to suggest that a properly formulated 
compressible formulation can be applied to flows that have 
been traditionally considered incompressible with essentially 
no penalty in accuracy or computational efficiency. It would 
seem that the computational barrier between the two flow 
regimes has been broken. The compressible flow formulation 
may be advantageous for low speed flows in which properties, 
including density, may vary significantly such as would occur 
in flows with heat transfer or chemical reactions. 

A second objective of the present paper is to examine the 
properties of a particular preconditioned formulation for some 
unsteady flows, namely the impulsively started driven cavity 
and the startup of pipe flow. Most of die early work on pre- 
conditioning concentrated on formulations for steady flows. 
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Tune-accurate preconditioning has been introduced relatively 
recently and the properties of such schemes for the Navier- 
Stokes equations have not been studied extensively over a 
range of flow configurations and Mach numbers. Whereas 
the merits of time-accurate formulations of Mach number 
preconditioned systems may have been demonstrated for the 
very low speed limit 7 , their accuracy and convergence have 
not been demonstrated over a very wide range of problems 
and Mach numbers. An additional minor goal of this paper 
is to supplement the discussions on preconditioning available 
in the literature with some physically based arguments and 
observations. 

In this study, the coupled strongly implicit procedure for 
solving the Navier-S tokes equations in primitive variables de- 
scribed in Ref. [8] was used as a starting point In the sections 
to follow, the modifications required to achieve convergence 
properties that are independent of Mach number at low speeds 
are first discussed in a one-dimensional context, then their 
implementation into the full two-dimensional Navier-Stokes 
equations is presented. Steady and unsteady results with and 
without preconditioning will be presented for the driven cav- 
ity problem. Results for the startup of pipe flow under a fixed 
total pressure will also be presented. The transient pipe flow 
problem was selected because it is one of the few physically 
attainable transient flows for which results are available for 
comparison. 


ANALYSIS 

Low Mach Number Limit 

It is well established that the incompressible form (including 
variable property versions) of the Navier-Stokes equations 
can be solved numerically in a straightforward manner. To 
reveal the source of the numerical difficulty that arises in 
solving the compressible form of the Navier-Stokes equations 
at low Mach numbers, we first examine the compressible 
Navier-Stokes equations in a form that should reduce to the 
incompressible but variable property form of the Navier- 
Stokes equations as the Mach number approaches zero. The 
density and other fluid properties may vary with temperature, 
but we anticipate that the density will become independent of 
pressure as the Mach number approaches zero. 

The key issues can be outlined sufficiently by using the 
one-dimensional Navier-Stokes equations. Nondimensional 
variables are defined as 
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where the variables with the tilde are the dimensional vari- 
ables, and the Mach number, M, is based on reference quanti- 
ties and the gas constant, M = Urcf /{jRTre } in order to 

recover the incompressible form of the equations in the most 


direct manner, primitive variables of u, p, and T will be used 
in the one-dimensional compressible formulation. There are 
no known disadvantages to the use of the primitive variables, 
particularly if the conservation law form of the equations is 
maintained. Substituting for density by using the ideal gas 
equation of state and utilizing primitive variables, p, u, and T, 
the conservation equations for mass, momentum and energy 
can be written 


where 
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The Reynolds and Prandtl numbers are defined as Re = 
PW «rc/£rc///*rc/, ft = C,ji/k. The Prandtl number of the 
fluid is assumed to be constant Now the consequences of the 
reference Mach number approaching zero will be considered. 
Note that p/RT = pjpr tJ and R = l/frM 2 ). We observe 
that the combination p/RT approaches a perfectly acceptable 
finite limit as M goes to zero. Note that E and E % will 
reduce to the traditional form of the incompressible Navier- 
Stokes equations as the Mach number approaches zero. The 
numerical solution will require linearization which can be 
accomplished by introducing Jacobian matrices which can be 
updated iteratively at each time step. This linearization is 
equivalent to solving 
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where the Jacobian matrices for the mass and momentum 
fluxes are given by 
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The linearization details for the viscous terms are being 
omitted since these terms do not influence the low Mach 
number issues under consideration here. In the above, R has 
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been replaced by 1 /{y M 2 ) in those terms that will vanish as 
M goes to zero. Notice that the variable property form of the 
equations usually applied for low speed flow of an ideal gas 
are recovered as M goes to zero. The mathematical nature of 
the time marching problem can be established by writing the 
system as 

S+wwg-wfc « 

and considering the eigenvalues of [At] -1 [A*]. The problem 
can be solved by a marching method if the eigenvalues are 
real. As M becomes small, [A t ] becomes ill-conditioned. 
That is, the determinant of [A,] becomes small (in fact, in the 
limit the first column and third row of [A t ] vanish) and errors 
are expected to arise in computing its inverse. In the limit as 
M goes to 0, the inverse of [A t ] does not exist and the system 
is singular. 

It is well-known that the eigenvalues of the matrix product 
[Aj -1 [A*] also provide an indication of the properties of the 
system. As the Mach number is decreased in the subsonic 
regime, the eigenvalues of the matrix product [At] -1 [-*4*] 
begin to differ more and more in magnitude (see for example, 
Turkel 1 ). As pointed out by previous investigators 2 ’ 3 the 
condition and degree of “stiffness” can be related to the 
relative magnitudes of the eigenvalues. When the eigenvalues 
differ greatly in magnitude, convergence to a steady state 
solution is usually slow, or for time-dependent solutions, the 
allowable time step becomes very small. This occurs because 
greatly varying signal speeds appear in the equations and the 
traditional solution schemes attempt to honor all of them, 
creating a “stiff” system. Since the Mach number has almost 
no influence on the physical characteristics of the flow at very 
low values of M (it effectively “cancels out” of the physics), 
it must be possible to devise a solution scheme in which the 
Mach number also “cancels out” over the range in which it is 
an unimportant physical parameter. 

It can be observed from [At] that as the Mach number 
approaches zero, the pressure is eliminated as an untaiown 
in the time derivative. This indicates that the solution to 
the incompressible equations cany no direct pressure history. 
Although the pressure field may not change much during a 
time step, in principle, the new pressure should depend solely 
on the new velocity field. Because the acoustic speed is 
effectively infinite in the incompressible limit, the pressure at 
a point may be influenced by the velocity field throughout the 
solution domain. A time-accurate compressible flow solution 
scheme must then provide a mechanism to permit the pressure 
field to be influenced in an elliptic fashion at each time step. 


where r is a pseudo-time and the conditioning matrix, [A p ] is 
given by 

[Ap] = Y RT 

i+MV^ 1 7M 4 p tl = 1 * ! 

Notice that [A p ] is formed from [A t ] by simply dividing 
the first column of [A t ] by 7 M 2 (equivalent to multiplying by 
the nondimensional gas constant R). Presumably, dividing by 
only M 2 would have the same effect. 

The hyperbolic system can be solved by advancing in 
pseudo-time until no changes occur at each physical time 
step. At that point, the time-accurate equations are satisfied. 
Obviously, this involves “subiterations”, but that is consistent 
with the observation that for a completely incompressible 
flow, the pressure field must be established at each physical 
time step with no direct dependence on a previous pressure 
field as was mentioned above. The hyperbolic system being 
solved is in r and x. It is similar to marching in pseudo- 
time to a “steady” solution at each physical time step. The 
addition of the pseudo-time term changes the eigenvalues of 
the hyperbolic system so that they are clustered closer together 
in magnitude, effectively at speeds closer to the convective 
speed. Eigenvalues for this system have been studied by 
Jorgenson 9 and found to be similar to those reported by 
Withington et al. 10 As the Mach number goes to zero, the 
ratio of the largest to smallest eigenvalues approaches 2.27 
for air with t=1.4. 

The second and minor source of difficulty with the coupled 
compressible formulation at low Mach numbers is related to 
the differences in the relative magnitudes of the nondimen- 
sional dependent variables. The magnitudes of the velocity 
and temperature remain of order one but the nondimensional 
pressure tends to increase without limit as the Mach number 
decreases. This permits round-off errors to become significant 
in the pressure gradient terms. Again, this is not a problem of 
physical origin. The present study indicated that no difficul- 
ties will usually be observed until the Mach number decreases 
to about 10 -5 or 10 -6 if double precision arithmetic is used. 
To completely eliminate the problem, a “gauge” or relative 
pressure 3 can be introduced so that differences in pressure 
(originating from pressure derivatives) become differences in 
the gauge pressure which can then be smaller in magnitude. 


Preconditioned Form for the Two-Dimensional 
Navier-Stokes Equations 


Overcoming the Problem at Low Mach Numbers 

To overcome the awkward (singular) mathematical situation 
that arises at low Mach numbers with the unsteady form of the 
coupled compressible equations and obtain a well-conditioned 
system, an appropriately formulated pseudo-time term can be 
added to the equations which vanishes at convergence at each 
physical time. Although the pseudo-time term can take many 
different forms and still be effective, simply adding a term 
of the same form as for the physical time term but with the 
Mach number removed from the terms causing the fatal ill- 
conditioning (the first column of the [A t ] matrix) is sufficient. 
The equations then become 


After replacing the density by pressure and temperature using 
the equation of state (p = p/RT ), the nondimensional form 
of the unsteady Navier-Stokes equations can be written in 
generalized nonorthogonal coordinates in a form applicable 
to both two-dimensional and axisymmetric flow as 

dQ(q] d(E(q) - E v (q}) 8(F(q) - F.(fl) _ „ 

at di dn 
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In the above, £ and t? are the spatial coordinates in the 
generalized system, and J is the transformation Jacobian. 
The parameter 6 is an index indicating the type of flow, two 
dimensional (6 = 0) or axisymmetric (S = 1). The parameter 
r is the nondimensional distance from the axis of symmetry 
in the axisymmetric cases. The nondimensional variables are 
defined as in the one-dimensional case with the addition of 
y = y/Lrtf and v = v/ur tf . 

The Navier-Stokes equations given above can also be 
written as 




where [A*], [A ( ], [A,], [£ £ ], [-B,] 
defined as 
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In particular, the Jacobian matrix [ A t ] is given by 
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To enable efficient numerical solution of the system at low 
Mach numbers, a term [A r ] §* is added to the equations where 
r is a pseudo-time and [A t ] is a matrix formed by dividing 
the first column of \A t ] by 7 M 2 (equivalent to multiplying 
the first column of [A t \ by the nondimensional gas constant 


R). The motivation for this was discussed previously. The 
equations to be solved numerically become 

r , W 

where the matrix [A T J is given by 


1 

T 

0 

0 

To* 

« 

T 

jSf 

0 

-rT* 

V 

T 

1 + iHLfcd 

0 

_e_ 

JL-. 
RT 
_E L_ 

1 1 

7 + 2 <vr 

RCpT 

RCpT 

2CpRT* 


Note that the choice of dependent variables does nothing to 
alter the conservation law form of the equations. 

The equations were discretized and solved numerically 
using the coupled strongly implicit procedure described by 
Chen and Pletcher.® The scheme utilized central differences, 
generally (some with deferred corrections®), and was first 
order in time. The equations were linearized by evaluating 
the Jacobian matrices at the previous pseudo-time step. Time 
accuracy was maintained by converging the system at each 
pseudo-time step, at which condition the pseudo-time term 
vanished. Thus, the converged equations at each physical 
time step satisfied the discretized Navier-Stokes equations. 
No smoothing was used for the cases shown here. 

For the driven cavity case, no slip boundary conditions 
were imposed on all boundaries. For the pipe flow cases, 
no slip conditions were used at the pipe wall. Symmetry 
conditions were imposed at the centerline. At outflow, the 
static pressure was specified and all other dependent variables 
extrapolated from the interior. 

In this study it was desired to simulate the startup of pipe 
flow as realistically as possible. In a typical application, 
a fixed total pressure is established and maintained at the 
entrance to the pipe and a valve at the discharge plane is 
opened. As the flow accelerates from rest, the static pressure 
of the moving fluid drops. The static pressure and velocities 
at the pipe inlet and throughout the length of the pipe vary 
with time until a steady flow is ultimately established. The 
value of the flow rate finally established will depend upon the 
difference between the total pressure fixed at the inlet and the 
fixed discharge pressure and the characteristics of the pipe. 

In the numerical simulations, the static pressure at the 
inlet was extrapolated from the interior and the velocities 
determined in order to maintain a constant total pressure. In 
addition, a constant total temperature was specified and the 
streamwise derivative of the normal component of velocity 
was taken as zero. That is, the conditions specified at the inlet 
are 
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where Pt, Tt, are specified and pi t j is extrapolated from 
the interior. The boundary conditions were implemented in 
an implicit fashion. The steady flow profile observed at the 
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Table 1: Driven Cavity, Re = 100, 19 x 19 Grid, Steady Row, 
Preconditioned ,A t = 10 3 , A t P — 0.3 


Mach No. 

10 _i 10" 3 10- 1 0.2 0.4 0.8 1.0 1.2 

No. of Iter. 

'52 52 52 52 52 52 53 53 


Table 2: Driven Cavity, Re=100, 19 x 19 Grid, Steady Row, 
No conditioning 


Mach No . 

.01 .1 .1 .1 .2 .2 .4 .4 .4 .8 

At 

.00025 .02 .025 .03 .05 .1 .15 .2 .25 .3 

No. of Iter. 

2123 355 302 264 182 138 82 68 97 51 


inlet of a two-dimensional channel for such a formulation 
corresponded well with the results reported by Wang and 
Longwell. 11 


RESULTS 


Impulsively Started Driven Cavity 
Steady State Solution 


Prior to undertaking detailed time dependent calculations, 
the behavior of the preconditioned system was evaluated and 
compared with unconditioned results over a range of Mach 
numbers. The computational effort in terms of iterations 
to achieve steady state is indicated for the preconditioned 
scheme in Table 1. Notice that the effort is essentially 
independent of Mach number and that the same pseudo-time 
step can be used over the entire range of Mach numbers 
considered. Many of these calculations were repeated on a 
finer 39x39 grid for Mach numbers of 0.4 and lower and the 
same trend was observed, namely, that the total number of 
iterations needed to achieve steady state was constant. For 
the finer grid, 100 iterations were required. The steady state 
solution was identical for Mach numbers of 0.1 or less. 

The computational effort required to obtain the steady flow 
solution without preconditioning is shown in Table 2. In these 
calculations it was not possible to use the same time step over 
a wide range of Mach numbers. Some experimentation was 
employed with the time step, but it was not easy to predict 
the optimum time step for a given Mach number. Note 
that the number of iterations was not reduced to the level 
achieved with the preconditioning until the Mach number 
reached 0.8. Convergence was determined for these steady 
calculations when the residual based on the norm of all 
variables 10 reached 10“ 5 . The steady flow solution details 
were virtually independent of Mach number for M less than 
about 0.1 The grid stretching was the same for all of the 
steady cases. 

Unsteady Solution 

For the unsteady driven cavity cases, the fluid was initially 
at rest and the pressure field was uniform. The unsteady 
motion started as the cavity lid was impulsively brought to a 


Table 3: Driven Cavity, Re = 100, 19x19 Grid, Unsteady 
Row, Preconditioned 


Mach No. 

lO- 1 

10~ 4 

10" a 

10" J 

10" 2 

10 _1 

0.2 

At 

0.05 

0.05 

0.05 

0.1 

0.05 

0.05 

0.05 

At p 

0.3 

0.3 

0.3 

0.3 

0.3 

0.3 

0.3 

Total Iter. 

651 

652 

652 

428 

652 

626 

538 

Sub-iter. 

4 

4 

4 

4 

4 

4 

4 


Table 4: Driven Cavity, 39x39 Grid, Unsteady Row, Pre- 
conditioned 


Mach No . 

.001 

.001 

.1 

.2 

.2 

.4 

.001 

.1 

.2 

.4 

Re No. 

100 

100 

100 

100 

100 

100 

400 

400 

400 

400 

A* 

.025 

.025 

.025 

.025 

.025 

.025 

.05 

.05 

.05 

.1 


.3 

.15 

.15 

.1 

.15 

.1 

.15 

.15 

.075 

.02 

Total Iter. 

967 

980 

926 

938 

881 

854 

1599 

1554 1697 

2232 

Sub-iter. 

4 

4 

3 

3 

4 

3 

4 

4 

4 

10 


constant velocity. The computational effort for the unsteady 
cases were generally independent of Mach number at the 
low end of the range considered when conditioning was 
employed. Table 3 indicates the total number of iterations 
for cases computed on the coarsest grid. The total iterations 
includes the subiterations required to achieve convergence to 
10 -5 at each time step. The average number of subiterations 
needed to achieve convergence at each pseudo-time step is 
also indicated in the table. 

On the coarsest grid, 19 x 19, the total number of iterations 
for the unsteady cases is seen to reach nearly a constant value 
of 652 as the Mach number decreases below 0.1. On the 
39x39 grid, more variation was observed in the total number 
of iterations as the Mach number was decreased while all other 
parameters remained fixed, although the variation was less 
than 10%. A possible explanation for this is as follows. The 
steady flow cases were computed using a very large physical 
time step. The result was that the physical time derivative 
term had no influence on the solution. Only the pseudo- 
time term influenced the numerics and the preconditioning 
effectively suppressed the influence of Mach number from 
that term. Thus, it is not surprising that the solutions become 
virtually independent of Mach number at the low end, i.e., 
as the Mach number approaches zero. For true transient 
calculations, however, the physical time term remains, and 
the time step needs to be sufficiently small to achieve time 
accuracy. The physical time term behaves like a source term 
during the pseudo-time iterations, and the size of the source 
term varies with Mach number. Thus, there may be no reason 
to expect the total iteration count to be rigidly independent 
of Mach number for time accurate calculations. This aspect 
of time-accurate preconditioned methods has not received 
much attention until now. Only a few unsteady cavity cases 
were computed without preconditioning. The computational 
effort required for the unconditioned scheme is summarized 
in Table 5. 

Note that the number of iterations increases significantly 
as the Mach number decreases below 0.2. Note also that 
subiterations were required at each time step even without 
the use of pseudo-time. The iterations were required to 
remove the linearization errors and converge the pressure 
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Table 5: Driven Cavity, Re = 100, Unsteady Flow, No 
Conditioning 


Mach No. 

0.01 

0.1 

0.2 

0.2 

0.4 

0.6 

Grid 

39x39 39x39 19x19 39x39 19x19 19x19 

At 

0.001 

0.01 

0.05 

0.025 

0.1 

0.05 

Total Iter. 

5288* 

1883 

905 

1167 

423 

472 

Sub-iter. 

4 

4 

5 

4 

4 

3 


•This case was only computed to t= 1.34 whereas steady state is achieved 
at approximately 6 time units. 
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Figure 1: Behavior of drag force at early times, impulsively 
started cavity. 
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Figure 2: Time history of predicted u component of velocity 
at cavity center, Re=100, impulsively started cavity. 



Figure 3: Time history of predicted u component of velocity 
at cavity center, Re=100, impulsively started cavity. 


field. The number of subiterations can be seen to decrease 
as the Mach number increases. At Mach numbers of 0.01 
and below, the transient results were in good agreement 
with those obtained by others 12 using purely incompressible 
formulations. Figure 1 compares the variation of the drag 
coefficient with time obtained with the present preconditioned 
compressible scheme with the computational results obtained 
by Soh and Goodrich. 12 The nondimensional time, f, in 
the figures for the cavity flow calculations is defined as 

Many aspects of steady driven cavity flow have been re- 
ported in the literature, particularly velocity profiles. Some 
aspects of unsteady incompressible driven cavity flow partic- 
ularly the time history of the flow patterns have been reported 
by Soh and Goodrich 12 . These will not be repeated here. 

The time history of the computed u component of velocity 
at the center of the cavity is given in Fig. 2 for M=0.001, 0.01 
and 0.1. The computations were carried out on a 39x39 grid. 
No conditioning was used for the results reported for M=0.0 1 . 
At first glance it is easy to conclude that the time history is 
essentially the same for these three Mach numbers. Closer 
inspection reveals a slight wiggle in the velocity for M=0. 1 at 
small nondimensional times. This was easy to ignore at first 
in this study, but it is now believed that this behavior is a real 
part of impulsively started compressible driven cavity flows. 
The oscillation in the velocity at the center becomes more 
pronounced as the Mach number increases, and as indicated 
in Fig. 3, where results for M=0.2 and 0.4 are shown, it cannot 
be ignored or considered an aberration. 

A very similar behavior can be observed in the solutions ob- 
tained without conditioning. That is, the oscillatory behavior 
appears not to be a product of the preconditioned formula- 
tion. This is indicated in Fig. 4 where both preconditioned 
results and results obtained without conditioning are shown 
for Re=100 on the 39x39 grid. The two results are nearly 
identical. Also shown in Fig. 4 are preconditioned results 
obtained on a 79 x79 grid. The oscillatory results are more 
sharply resolved on the finer grid, but the qualitative picture 
remains the same. 

A possible explanation for this behavior is as follows. The 
impulsive start of the cavity lid creates a singularity in the flow 
field. That is, an absolutely impulsive start is not physically 
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Figure 4: Effect of grid refinement on time history of u 
component of velocity at cavity center, Re=100, impulsively 
started cavity. 
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possible. Any real start of the motion must occur over a 
finite interval of time. In their incompressible work, Soh and 
Goodrich pointed out a square root singularity in the drag on 
the moving lid. The fact that the drag curve in Fig. 1 islinearat 
small times tends to verify this singular behavior. To resolve 
the flow accurately at very small times, finer and finer spatial 
resolution is required. For purely incompressible flow, the 
effects of the singularity are transmitted only by convection 
and diffusion. The initial pressure has no influence on the 
solution. Acoustic signals are transmitted instantaneously so 
any influence that the abrupt start has on the pressure field is 
mitigated by the rapid propagation of signals. 

The magnitude of the reference Mach number in this prob- 
lem is precisely the nondimensional time required for an 
acoustic wave to move the length of the cavity if propagating 
at the sonic speed corresponding to the reference temperature. 
For calculations at very low Mach numbers, an acoustic wave 
could make several traverses across the cavity during the first 
time step. This leads to the results seen in Fig. 2 which is 
exactly the behavior of solutions to the incompressible equa- 
tions. At very low Mach numbers, the physical time term 
(source term for solution advancement in pseudo-time) is very 
small, and the effect of pressure history is correspondingly 
small. For the calculations at M=0.1, the time required for 
acoustic waves to traverse the cavity is about 0.1 nondimen- 
sional time units. This time interval is now larger than the 
computational time step required to maintain temporal accu- 
racy and a slight oscillation can be observed in Fig. 2. The 
oscillations can be seem more clearly in Fig. 3. The period of 
the initial oscillations observed in Fig. 2 is approximately 0.4 
for the M=0.2 results and 0.8 for the M=0.4 results. The pe- 
riods are approximately the time required for acoustic waves 
to propagate across the cavity and back. The oscillations are 
seen to damp out as the solution tends toward steady state. 
The pressure at the center of the cavity is observed to possess 
a similar oscillatory pattern as can be seen in Fig. 5. The p re f 
in Fig. 5 is the pressure value at the center of the bottom plate 
of the cavity. 

The oscillatory pattern is seen to sharpen somewhat when 


Figure 5: Time history of predicted pressure at cavity center, 
Re=100, impulsively started cavity. 



Time, t 

Figure 6: Effect of time step on the time history of u 
component of velocity at cavity center, Re=100, impulsively 
started cavity. 
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Figure 7: Time history of u component of velocity at cavity 
center, Re=400, impulsiveUy started cavity. 

a smaller physical time step is employed as can be seen in 
Fig. 6. This is reasonable since the phenomenon is believed 
to originate with the impulse associated with the start of the 
lid. The smallest time step that will permit the viscous region 
to reach the first grid point from the wall (for the 39 x 39 grid) 
in an analogy with the impulsively started flat plate is about 
0.0125 corresponding to the solid line results in Fig. 6. The 
trends observed for a Reynolds number of 400 are somewhat 
similar in that oscillatory behavior was observed for Mach 
numbers of about 0.2 or larger. The flow at Re=400 is actually 
quite different and the time history is quite interesting 12 . For 
example, the history of the velocity at the center is not 
monotonic and the nondimensional time required for the flow 
to reach a steady state is 2-3 times as long as for the flow 
at Re=100. Figure 7 shows the history of the x-component 
velocity for calculations at Mach numbers of 0.001, 0.2, and 
0.4. Oscillations are detectable in the M=0.2 and 0.4 results. 
Again, the initial periods of the oscillations are about 0.4 and 
0.8, respectively, although this is not evident from the plots. 


Startup of Pipe Flow 

A goal of the present study has been to evaluate the time 
accuracy of preconditioned schemes with a view toward 
determining the merits of using a single scheme for simulating 
both incompressible and compressible phenomena. The 
startup of flow in a pipe is a real transient phenomena that 
can be observed in the laboratory with greater ease than the 
impulsively started cavity considered above. Pipe flow is a 
pressure driven flow whereas the cavity is a shear driven flow 
and it was speculated that compressibilty effects might play a 
more dominant role in the physics of transient pressure driven 
flows. 

Incompressible starting pipe flow has been studied com- 
putationally through solutions to theNavier-Stokes equations 


by Anderson and Kristoffersen 13 and through solutions to the 
boundary-layer equations by Ramin and Pletcher. 14 For an 
infinitely long pipe (no entrance effects) and incompressible 
flow, the governing equations for the startup problem be- 
come linear and a solution in the form of Bessel functions 
was proposed by Szymanski. 15 For pipes of finite length, 
a flow development region exists and the flow is governed 
by the nonlinear Navier-Stokes equations. Except for small 
Reynolds numbers, the boundary-layer approximation can 
provide a model that provides reasonably accurate solutions 
over most of the development length under incompressible 
assumptions. The full Navier-Stokes equations were used in 
the present study. 

In the transient pipe flow problems considered here, the 
total pressure was considered to be fixed at the pipe inlet. The 
static pressure at the pipe discharge was also fixed and this 
pressure difference [H=P to t«i (at inlet)- P,tatic (at discharge)] 
became part of the problem specification. The appropriate 
characteristic velocity for this transient flow is not obvious. 
The final steady flow velocity (or mass flow rate) is not known 
a priori. One physically relevant characteristic velocity is 
given by U p =(H/p )*/ 2 , the maximum velocity that could 
be achieved from an application of Bernoulli’s equation on 
a lumped (stream tube) basis to this flow. A second one 
is the average velocity that would exist in incompressible 
laminar fully developed flow consuming the pressure drop H, 
V ,=££i 

>n J 32ftL ' 

Otis 16 , working with incompressible flow assump- 
tions, proposed the use of the nondimensional parameter 
E=^ ( y ) 2 , where Re p is the Reynolds number based on 
diameter and U p , (Otis and others used M for this parameter 
but M is not appropriate here) to characterize this startup 
flow. Note that E=0 corresponds to flow in an infinitely 
long pipe. The parameter E can be thought of as the ratio 
of the time scales for cross stream diffusion and streamwise 
convection. The value of E will also determine how closely 
the steady state flow rate will approach the flow that would be 
predicted by considering H to be the static pressure drop for 
fully developed laminar flow. The larger the value of E, the 
smaller the flowrate compared to the flow that would result 
if the H were consumed by pressure drop in fully developed 
laminar flow. In other words, the larger the value of E, the 
greater the entrance effect and its corresponding momentum 
pressure drop. Nearly all of the previous literature on this 
subject dealt with incompressible flow. All computations 
were made using air as the working fluid, isothermal wall 
conditions with the Prandtl number set equal to 1.0. The 
Mach number reported for the pipe results was based on the 
Vinf reference velocity defined above and the sonic velocity 
based on the inflow stagnation temperature. 

It was of interest first to determine if the present com- 
pressible formulation could duplicate all the expected physics 
for a purely incompressible flow. To this end, computations 
were made for the startup problem using a Mach number of 
1.386 xl0~ 3 . Calculations were made for E=0.084, 0.108, 
and 0.5. Two grids were employed for the E=0.108 case. 
To simulate incompressible flow calculations correctly, at all 
time steps the mass flow rate should be everywhere the same 
in the pipe since the density of an incompressible fluid would 
not change with time. This was achieved to within less than 
0.5%. In terms of the nondimensional equations, the small 
value of Mach number is effectively suppressing the physical 
time term in the continuity equation. 

Figure 8 shows the time history of the mass flow rate for 
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Figure 10: Predicted pipe discharge velocity distributions at 
various times, M=0.001386,E=0.108. 


Time, t* 

Figure 8: Time history of predicted (M=0.001386) pipe 
discharge rate for several values of E. 



Figure 9: Predicted pipe inlet velocity distributions at various 
times, M=0.001386, E=0.108. 


the M=1.386xl0 -3 cases computed. The nondimensional 
time, <*, in the figures for the pipe flow calculations is 

defined as t*=t/(^). The flowrate in the figure has been 
nondimensionalized with the flowrate that would be achieved 
for a fully developed incompressible laminar flow consuming 
the same pressure drop in the same tube, i.e., based on an 
average velocity of V in} defined above. The present results 
are seen to be in good agreement with the numerical results of 
Anderson and Kristoffersen 13 , especially for E=0.108. Grid 
refinement for that case made little difference in the results as 
can be seen in Fig. 8. The flow at E=0.108 reaches the fully 
developed state with the centerline velocity equal to twice the 
average velocity. For the E=0.5 case, the centerline velocity 
was 3% short of the fully developed value. Velocity profiles 
at inlet and at the discharge are shown in Figs. 9 and 10. It 
is of interest to note that the inlet velocity distribution is not 
uniform under the present method of fixing the total pressure 
at inlet. 

The transient behavior of flow was strikingly different for 
the higher Mach number cases computed. At early times, the 
fluid in the pipe “decompressed” allowing flow to discharge 
from the pipe while the inlet flow was zero. This can be 
seen in Fig. 11. At the earliest time shown, one third to 
one half of the fluid is at rest as fluid is forced out the 
discharge at flow rates sometimes higher than the steady state 
values. The variation of the discharge rate with time for 
three cases, M=0.1386, M=0.2772 and M=0.3224 are shown 
in Fig. 12. This is in contrast to the results of the lowest 
Mach number case (effectively incompressible) shown in Fig. 
8. Even though the steady flow behavior of all of these cases 
may be close to that attributed to incompressible flow, the 
transient behavior is greatly different. A design based on 
incompressible results shown in Fig. 8 may not work for 
flows at speeds as low as approximately Mach 0.1. 

The time history of the centerline pressure is shown in 
Fig. 13 for the two highest Mxh number cases. The pressure 
remains at nearly the total pressure for the majority of the pipe 
for the earliest time shown. In contrast, for an incompressible 
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Figure 11: Predicted distribution of flow rate along pipe at 
early times, E=0. 108. 



Figure 12: Predicted time history of pipe discharge at three 
values of Mach number, E=0.108. 



Figure 13: Predicted distribution of centerline pressure along 
pipe at early times, E=0.108. 

case (not shown in Fig. 13) the earliest pressure would be 
more like the large time curves shown in Fig. 13. For the pipe 
flow cases, the nondimensional time required for an acoustic 
wave, moving at the sonic speed corresponding to the inlet 
total temperature, to travel the length of the pipe is given by 
M/( 16 E). 

The computing effort required for the pipe cases varied 
noticeably over the Mach number range. This is in contrast 
to the driven cavity case where the physics did not change 
much as the Mach number varied. In the pipe flow case, the 
flow was driven by pressure differences. At the lowest Mach 
number considered, essentially no motion could occur until 
the pressure field was established in line with the boundary 
conditions and flow occurred simultaneously throughout the 
entire computational domain. On the other hand, essentially 
no action occurred in a large portion of the pipe at early times 
for the higher Mach number cases. Constant time steps of 
t* =0.0133 were used for all results shown. Approximately 85 
time steps were needed to achieve the steady state conditions. 
Most of the calculations were carried out on a 30 x 10 grid but 
two calculations were made on 50x20 grids for evaluation of 
the effects of grid refinement. The convergence criteria for 
each time step was set at 10 -6 for the pipe flow cases. For 
steady pipe flow calculations, starting with the fluid at rest, 
718 iterations were required for the M=0.001386 case and 
933 for the M=0.2772 case both on a 50x20 grid. For the 
finest grid unsteady flow calculations, however, an average 
of 310 iterations was required for convergence at each time 
step for the M=0.001386 case whereas only an average of 61 
iterations was required per time step for the M=0.2772 case. 
This difference was not anticipated in light of the experience 
reported for the driven cavity transient calculations. Whether 
this is reasonable in light of the the complexity of the physics 
is presently under study. 


10 





CONCLUSION 

Generally, the coupled, preconditioned compressible for- 
mulation was observed to perform well over a range of Mach 
numbers for two transient test cases. For the driven cav- 
ity flow, the computational effort was nearly independent 
of Mach number especially at low values although this was 
more precisely true for the steady cavity calculations than 
for the transient ones. The impulsive start of the cavity 
lid was observed to promote an oscillatory behavior in the 
velocities and pressure that was not observed at the lowest 
Mach numbers. The preconditioned results, including the 
oscillatory velocities, were in good agreement with results 
obtained without preconditioning. 

The calculations of the startup of pipe flow obtained for a 
very small Mach number, M=0.001386, were in agreement 
with results obtained by Anderson and Kristoffersen 13 who 
solved the incompressible equations. For the higher Mach 
numbers considered, the transient behavior was entirely dif- 
ferent in that the flow started by decompression while the 
inlet mass flow remained essentially zero. That this should 
be true can be reasoned taking into account the time required 
for acoustic waves to penetrate from the discharge plane to 
the inlet and the compressibility of the fluid. Conservation of 
mass can be maintained without inflow as the density changes 
with time. The computational effort required to solve the 
lowest Mach number case for the pipe startup was substan- 
tially greater than that required for the higher Mach number 
cases. Such a dependency has not been observed for previous 
steady flow calculations. This issue is presently under study. 
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